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In this letter we present a flat histogram algorithm based on the pruned and enriched Rosen- 
bluth method (PERM). This algorithm incorporates in a straightforward manner microcanonical 
reweighting techniques, leading to "flat histogram" sampling in the chosen parameter space. As an 
additional benefit, our algorithm is completely parameter free and, hence, easy to implement. We 
apply this algorithm to interacting self-avoiding walks (IS AW), the generic lattice model of polymer 
collapse. 



Recently, there has been revived interest in fiat his- 
togram algorithms pj , which strive to evenly sample con- 
figuration space with respect to a chosen parameteriza- 
tion, e.g. microcanonical energy. These algorithms are 
particular implementations of "umbrella sampling" , in 
which the configuration space is sampled according to a 
given probability distribution, the so-called "umbrella". 
This umbrella distribution is generally chosen such that 
the whole configuration space of interest is accessible in 
one simulation. One major difficulty hereby is to find a 
suitable umbrella distribution. 

There has also been an exciting development in 
stochastic growth algorithms, which are based on the 
Rosenbluth and Rosenbluth algorithm 0. If this algo- 
rithm, which kinetically grows configurations, gets en- 
hanced by cleverly chosen enrichment and pruning steps 

0, one obtains the pruned and enriched Rosenbluth 
method (PERM), a powerful algorithm for, e.g., simu- 
lation of the polymer collapse transition. 

We present in this letter a new algorithm, flatPERM, 
which is a combination of these two types of algorithms, 

1. e. a flat histogram version of the pruned and enriched 
Rosenbluth method. 

As opposed to earlier work in this direction |5(, in 
which an iterative scheme similar to the multicanonical 
algorithm Q was used, we utilize the self-tuning capa- 
bilities of PERM directly. This leads to a considerable 
simplification of the algorithm. 

While flatPERM includes umbrella sampling ideas, it 
is stricly speaking not a multicanonical method, as multi- 
canonical sampling conventionally describes a particular 
iterative version of adaptive umbrella sampling, in which 
first an umbrella distribution is obtained iteratively and 
then a final simulation is performed with a fixed umbrella 
distribution. 

This letter is structured as follows: we first give a ped- 
agogical introduction to PERM, which will then allow us 
to introduce flatPERM as a seemingly trivial extension. 
As an application, we present simulations of interacting 
self-avoiding walks (ISAW) on the square lattice and the 
simple cubic lattice. We conclude with a description of 
further applications. 

We will consider a rather abstract setting of configura- 
tions with a certain size n, which is parameterized with 
an additional variable m. In general, one can even con- 



sider a set of variables rrij, but for pedagogical reasons 
we shall in this paper restrict ourselves to the case of m 
corresponding to an energy E = em. Both n and m are 
assumed to have non-negative integer values. Moreover, 
we will need the notion of "atmosphere" of a configura- 
tion, which is the number a of different ways to continue 
to grow this configuration, and is also a non-negative in- 
teger. 

While it is useful to present the algorithm in such an 
abstract setting, it may help the reader to keep the ap- 
plication to interacting self-avoiding walks (ISAW) on a 
regular lattice in mind. In this case, the size of the config- 
uration is the number of steps of the walk, the energy is 
the number of non-consecutive nearest-neighbor bonds, 
and the atmosphere is the number of non-occupied sites 
around the endpoint of the walk. If all the sites around 
the endpoint are occupied, then the atmosphere is zero 
and the walk cannot be continued. 

The basis of the algorithm is the Rosenbluth and 
Rosenbluth algorithm, a stochastic growth algorithm 
in which the configurations of interest are grown from 
scratch. The growth is kinetic, which is to say that 
each growth step is selected at random from all possi- 
ble growth steps. Thus, if there are a possible ways of 
growth, one selects one of them with probability p = 1 /a. 
As this number generally changes during the growth pro- 
cess, different configurations are thus generated with dif- 
ferent probabilities, and one needs to resort to reweight- 
ing techniques. 

It is advantageous to view this algorithm as approxi- 
mate counting, in which case the weights of the configura- 
tions serve as estimates of the number of configurations. 
To understand this point of view, imagine that we were 
to perform a complete enumeration of the configuration 
space. Doing this requires that at each growth step we 
would have to choose all the possible continuations and 
count them each with equal weight. If we now select fewer 
configurations, we have to change the weight of these con- 
figurations accordingly, in order to correct for missing out 
on some configurations. Thus, if we choose one growth 
step out of a possible ones, we replace a configurations 
with equal weight by one "representative" configuration 
with a- fold weight. In this way, the weight of each grown 
configuration is a direct estimate of the total number of 
configurations. 
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Let the atmosphere a n = a(ip n ) be the number of dis- 
tinct ways in which a configuration tp n of size n can be 
extended. Then, the weight associated with a configu- 
ration of size n is the product of all the atmospheres a k 
encountered while growing this configuration, i.e. 

n-l 

W=l[a k . (1) 

k=0 

After having started S growth chains, an estimator C^ st 
for the total number of configurations C n (the "infinite- 
temperature" partition function, in which all configura- 
tions appear with equal weight 1) is given by the mean 
over all generated samples ipn of size n with respective 

(i) 

weights Wn , i.e. 

C? = (W) n = ±J2 W n ] ■ (2) 

i 

Here, the mean is taken with respect to the total number 
of growth chains S, and not the number of configurations 
which actually reach size n. Configurations which get 
trapped before they reach size n appear in this sum with 
weight zero. 

The problem with the Rosenbluth and Rosenbluth al- 
gorithm is two-fold. Firstly, the weights can vary widely 
in magnitude, so that the mean can get to be dominated 
by very few samples with very large weight. Secondly, 
regularly occurring trapping events, i.e. generation of 
configurations with zero atmosphere, can lead to expo- 
nential "attrition" , i.e. exponentially strong suppression 
of configurations of large sizes. 

To overcome both of these problems, enrichment and 
pruning steps have been added to this algorithm, leading 
to the pruned and enriched Rosenbluth method (PERM) 
p|. The basic idea is that one wishes to suppress fluc- 

tuations in the weights Wn , as these should on average 
be equal to C n . Therefore, if the weight is too large, 
one enriches, i.e. one makes copies of the configuration 
and reduces the weight accordingly. On the other hand, 
if the weight is too small, one prunes probabilistically, 
and, in case the pruning is unsuccessful, continues grow- 
ing with appropriately increased weight. Note that 5, 
and therefore the expression @ for the estimate C^ st , is 
not changed by either enriching or pruning steps. 

We need to specify enrichment and pruning criteria 
as well as the actual enrichment and pruning processes. 
While the idea of PERM itself is straightforward, there is 
now a lot of freedom in the precise choice of the pruning 
and the enrichment steps. The "art" of making PERM 
perform efficiently is based to a large extent on a suitable 
choice of these steps. Distilling our own experience with 
PERM, we present here what we believe to be an efficient, 
and, most importantly, parameter free version. 

In contrast to other expositions of PERM (e.g. .7)), we 
propose to prune or enrich constantly, to enable larger 
exploration of the configuration space (the motivation 
will become clear once flatPERM is introduced below). 



Define r as the ratio of weight and estimated number of 
configurations, r = Wn /C® st . Then we enrich if r > 1 
and prune if r < 1. Moreover, the actual pruning and 
enrichment steps are chosen such that the weights are set 
as closely as possible to C^ st to minimize fluctuations: 

• r > 1 — > enrichment step: 

make c = min([? , J,a n ) distinct copies, each with 

weight \W^ ] (as in nPERM @). 

• r < 1 — > pruning step: 

continue growing with probability r and weight 
C® s * (i.e. prune with probability 1 — r). 

Note that we perform pruning and enrichment after the 
configuration has been included in the calculation of C^ st 
and is used for weights during the next growth step. Ini- 
tially, the estimates C^ 3t can of course be grossly wrong, 
as the algorithm knows nothing about the system it is 
simulating. However, even if initially "wrong" estimates 
are used for pruning and enrichment, in all applications 
considered the algorithm can be seen to converge to the 
true values. It is, in a sense, self-tuning. 

At this point it is now straight-forward to change 
PERM to a thermal ensemble with inverse tempera- 
ture — 1/ksT and energy E = em by multiply- 
ing the weight with the appropriate Boltzmann factor 
exp(— (3E), which leads to an estimate of the partition 
function Z n {j3) of the form 

Z° st ((3) = (Wexp(-{3E)) n . (3) 

The pruning and enrichment procedures are changed ac- 
cordingly, replacing W by W exp(—(3E) and C^ st by 

Z* st ((3), and using r = W$ exp(-(3E$)/Z* st ((3). (It 
is in this setting that PERM is usually described.) 

Alternatively, however, it is also possible to consider 
microcanonical estimators for the total number C n . m of 
configurations of size n with energy m (i.e. the "density 
of states" ) . An appropriate estimator C^ s ^ is then given 

' (i) 

by the mean over all generated samples ipn'm of size n 

(i) 

and energy m with respective weights Wn.m, i.e. 

C e n %=(W) n , m = ^J2 W ^k- ( 4 ) 

i 

Again, the mean is taken with respect to the total number 
of growth chains S, and not the number of configurations 
S n ,m which actually reach a configuration of size n and 
energy m. The pruning and enrichment procedures are 
also changed accordingly, replacing C n by C n ^ m and using 

r = Wn t m/Cn,m- 

At this point it is worth pointing out that the pruning 
and enrichment criterion for PERM leads to a roughly 
constant number of samples being generated at each size 
n for PERM. In fact, one can motivate the given pruning 
and enrichment criteria by stipulating that one wishes 
to have a roughly constant number of samples, as this 
leads to the algorithm performing an unbiased random 
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walk in the configuration size. Correspondingly, in the 
flat-energy version the algorithm performs an unbiased 
random walk in both size and energy of the configuration, 
and we obtain a roughly constant number of samples for 
each size n and energy to. 

It is because of the fact that the number of samples is 
roughly constant in each histogram entry, that this algo- 
rithm can be seen as a "flat histogram" algorithm, which 
we consequently call flatPERM. In hindsight in becomes 
clear that PERM itself can be seen as a flat histogram 
algorithm, at it creates a roughly flat histogram in size 
n. Recognizing this led us to the formulation of this al- 
gorithm in the first place. 

At this point we return to our discussion of the pruning 
and enrichment strategies. For PERM it may be more 
advantageous to allow for a high diffusivity along the size 
n. This is done by minimizing pruning and enrichment, 
i.e. at the cost of allowing larger weight fluctuations. For 
flatPERM on the other hand we need to achieve diffusive 
behavior also with respect to the energy variable to. Pre- 
cisely this is achieved by allowing for much enrichment 
(and thus necessarily also pruning) , as each set of config- 
urations enriched at size n will contribute to a range of 
different energies to at size n + 1. 

Even though we view flatPERM as a flat histogram al- 
gorithm, we have not yet explicitly conditioned the prun- 
ing and enrichment with respect to the local number of 
samples S n>m created at each size n and energy to. This 
can be done by multiplying r by S/ S ntm , as this enhances 
enrichment if S n ^ m < S, i.e. when there are too few sam- 
ples, and vice versa. 

One general problem with PERM is that enrichment 
generates large correlation between samples, so that it 
would be useful to replace the number S n m of samples 
by a number S%frn °f effectively independent samples, 
thereby taking account of some autocorrelation time. 
This can be done hcuristically by considering the average 
number of independent growth steps in a configuration. 
If the last enrichment has occurred at then a 

configuration of size n and energy m has riind = n — n enr 
independent steps. The frequency of independent steps 



in a configuration is w 



ni n d/n. We thus use 
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(5) 



for the pruning and enrichment criterion in our algo- 
rithm. 

Now that one has a flat histogram method of PERM 
for the whole range of energies, one can easily modify 
the algorithm further to sample only a selected range by 
dividing r by a "profile shape" f n , m , leading to pruning 
whenever /„.,„ is close to zero. This is advantageous if 
one only wants to explore a restricted region of parameter 
space. 

Even though the algorithm described thus far is free of 
parameters, there is a technical problem due to the fact 
that initial weights are grossly wrong, which can lead 



to overflow problems. This can easily be overcome by 
initially restricting the maximal size of grown configura- 
tions, e.g. by limiting the maximal size by the number 
of started growth chains, 



n < cS 



(6) 



The relevant number of growth chains at size n is thus 
reduced to S — n/c. For IS AW we find that a choice of 
delay of c ss 10 is sufficient. 
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FIG. 1: Logarithm of the number of configurations C n ,m 
versus internal energy m/n and length n for IS AW on the 
square lattice (above) and simple cubic lattice (below). 

Averages of observables Q defined on the set of config- 
urations can now be obtained by storing weighted sums 
of these observables, from which one obtains 
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(QW) n , m _ Qn!mWj% 

(W) n>m 



(7) 



These can then be used for subsequent evaluations. For 
instance, the expectation value of Q in the canonical en- 
semble at a given temperature /3 can now be computed 
via 



£ m QiM£e*P(-/^ 



(8) 



We have implemented this algorithm for interacting 
self-avoiding walks on the square and simple cubic lat- 
tice. In both two and three dimensions, we have sim- 
ulated walks up to length n = 1024. Figures ^ an( i [21 
clearly show the strength of the method. Figure ^ shows 
the number of configurations C rhm , which vary over sev- 
eral hundred orders of magnitude. This range would have 
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FIG. 2: Number of generated samples versus internal energy 
m/n and length n for IS AW on square lattice. The top graph 
shows the number of samples S n ,m for a simulation in which a 
flat histogram for S n , m was created, whereas the middle and 
bottom graphs show the number of samples S n ,m and the ef- 
fective number of samples S„{m, respectively, for a simulation 
in which a flat histogram for S^m was created. 

been inaccessible during one simulation with any canon- 
ical method. As an aside, we note that we had to rescale 
weights during the run to avoid overflow Addition- 
ally, we chose a delay of c = 10 to stabilize the algorithm. 

As can clearly be seen in Figure |3 the number of 
samples generated is roughly constant, irrespective of 
whether one flattens with respect to S n ,m (top graph) or 
S%fm (middle and bottom graphs). Only results for the 
square lattice are shown, as we find the same behavior for 
the simple cubic lattice. Using S%fm f° r the flatness cri- 
terion leads moreover to an increased sampling of walks 
with very few and very many interactions, thus overcom- 
ing the usual difficulty of obtaining configurations with 
a large number of interactions. (The largest energy state 
gets repeatedly sampled in simulations in both dimen- 
sions.) 

Once the simulations have been completed, ther- 
modynamic quantities of interest can easily be com- 
puted. As an example, the specific heat curves C n = 
kB(f3e) 2 a 2 (m)/n near the transition temperature are 



FIG. 3: Normalized fluctuations a 2 (m)/n versus inverse tem- 
perature [3 = l/ksT for IS AW on the square lattice (above) 
and the simple cubic lattice (below) at lengths 64, 128, 256, 
512, and 1024. The curves for larger lengths are more highly 
peaked. The vertical lines denote the expected transition tem- 
perature at infinite length. 



shown for both systems in Figure 

We have also tested our algorithm on the HP model 
[To), which is a toy model for proteins. It consists of 
a self-avoiding walk with two types of monomers along 
the sites visited by the walk, which are either hydropho- 
bic (type H) or polar (type P). One considers monomer- 
specific interactions, mimicking the interaction with a po- 
lar solvent such as water. The interaction strengths are 
thus chosen such that HH-contacts are favored, e.g. by 
choosing ehh = — 1 and chp = tpp = 0. The central 
question is to determine the density of states (and to find 
the ground state with lowest energy) for a given sequence 
of monomers. For various sequences taken from the lit- 
erature we have confirmed previous density of states cal- 
culations and obtained identical ground state energies. 
The sequences we considered had n = 58 steps (3 di- 
mensions, E m i n = —44) and n = 85 steps (2 dimen- 
sions, E min = —53) from 0, and n = 80 steps (3 di- 
mensions, E m i n = —98) from [TT|. We studied also a 
particularly difficult sequence with n — 88 steps (3 di- 
mensions, E m i n = —73) from 0, but the lowest energy 
we obtained was E = —69. However, in this case no 
lower energy has been found with any other PERM im- 
plementation either, see 0. 

To summarize, we have presented a new algorithm, 
flatPERM, which is a flat histogram version of a stochas- 
tic growth algorithm. This algorithm can in principle be 
applied to any statistical mechanical system for which 
configurations can be grown in a well-defined manner. 
Next to the presented applications of linear polymers, 
this algorithm can be applied to more complicated sys- 
tems, such as lattice models of branched polymers |13| . 
Extensions to models with two energy parameters are in 
preparation, e.g. to the problem of adsorbing interact- 
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